clear;clc;

FolderList      =   {'Shared','Temp', ...
                     'SteadyState','Equations','Graph'};
for ii=1:length(FolderList)
    addpath(genpath(FolderList{ii}));
end

%% Setup
PP          =   Setup_PP(struct('KAPPA',0.06/4));

%% Steady State
PSI_0       =   1;
DiscFactor_0=   0.02;
SS          =   SteadyState_Solver(PP,PSI_0,DiscFactor_0);


%% Construct the Model Structure
MODEL       =   Setup_MODEL(PP,SS);

% Check the Steady State in Equation Blocks
Res_NK      =   Equ_NK(PP,SS,0);
Res_HH      =   Equ_HH(PP,SS,MODEL,0);

if max(abs([Res_NK;Res_HH]))<1e-5
    fprintf('Steady State is Correct!\n');
else
    error('Steaty State is not Correct, check it  ...');
end


%% Compute the Jacobian
EquJac      =   struct();
EquJac.HH   =   Equ_HH(PP,SS,MODEL,1,[],[],1);

%% IRF

% Flexible Exchange Rate
PP.Taylor_Pi    =   1.5;
EquJac.NK   =   Equ_NK(PP,SS,1);
fir_ord     =   JacAssemble(MODEL,EquJac);
tic;
[gx,hx,gxhx_ExitFlag,gxhx_Info]=gx_hx(fir_ord,1);
toc;
if gxhx_ExitFlag~=1
    warning('Perturbation Solution is problematic');
    return
else
    solution    =   struct('gx',gx,'hx',hx);
end
solution    =   struct('gx',gx,'hx',hx);
StdVec      =   [1 -1/4]'*0.01;
IRF         =   IRF_1order(MODEL,solution,StdVec);


